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Abstract. We introduce efficient numerical methods for generic HJM equa- 
tions of interest rate theory by means of high-order weak approximation schemes. 
These schemes allow for QMC implementations due to the relatively low di- 
mensional integration space. The complexity of the resulting algorithm is 
considerably lower than the complexity of multi-level MC algorithms as long 
as the optimal order of QMC-convergence is guaranteed. In order to make 
the methods applicable to real world problems, we introduce and use the set- 
ting of weighted function spaces, such that unbounded payoffs and unbounded 
characteristics of the equations in question are still allowed. We also provide 
an implementation, where we efficiently calibrate an HJM equation to caplet 
data. 

1. Introduction 

The Heath- Jarrow-Morton equation (HJM-equation) of interest rate theory f[26j: 
see [2J [5J [T3] for expositions) is a stochastic partial differential equation (SPDE) 
on the state space of forward rate curves, which is flexible enough to describe com- 
plicated dynamical features such as non-constant (local or stochastic) volatility, 
non-constant correlation, or jumps, or dependence structures. An analysis of geo- 
metric properties was performed in [16] . As forward rate curves already encode all 
the market's information on default-free bond prices, it only remains to estimate 
volatilities either from the time series or from option prices or from both of them. 
For this purpose it is required that the numerical treatment of the HJM-equation 
can be performed efficiently: it is the purpose of this article to actually show that 
efficient numerical methods for the HJM-equation are at hand, how to construct 
and how to implement them. 

In the case of generic SPDEs we usually neither have sufficient analytical infor- 
mation on the marginal's distribution, nor on its Fourier-Laplace transform, nor 
its short-time asymptotics. We are therefore forced to apply simulation techniques 
to approximate the random variables in question and we face two main sources of 
problems in such a procedure: 

1.1. Discretization error. Numerical weak or strong approximation schemes with 
probabilistic flavor are built upon stochastic Taylor expansion and its iteration along 
n steps due to the Markov property. Depending on the local error of the method 
this leads (at least for some class of test functions) to a global error 0(l/n s ), which 
is called error of order s. The method is called high order method if s > 1, and 
standard or low order method otherwise. There are schemes, e.g. cubature methods 
pTTl 1521 1551 or splitting methods [55J EZ] ; that substantially increase s and therefore 
reduce the global discretization error for a fixed number of discretization steps n. 
When applying the theory of weighted spaces we can also enlarge the sets of test 
functions and generic equations, which can be treated by the discretization method. 

1.2. Integration error. Having discretized the SPDE problem we still have to 
evaluate the random variables involved in each local step, which usually leads to a 
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numerical integration problem on some M. dn , where d is the fixed number of dimen- 
sions which are needed for each local discretization step. Here we can apply three 
approaches: (deterministic) numerical integration, Monte-Carlo algorithms (MC) 
or Quasi-Monte-Carlo algorithms (QMC). Due to the n-dependence of integration 
space we do not try a direct numerical integration method, even though we have 
some hope that such an approach could possibly work. MC algorithms lead to 
integration errors Q(1/\K), where K denotes the number of integration points, 
whereas QMC algorithms lead to integration errors approximately 0(1/ K). In both 
cases the integration error dominates the total error asymptotically, which can be 
seen by complexity analysis. Let us be more precise on this: we assume that 0(n) 
operations are performed to calculate the value of the functional which we intend 
to integrate. Here we tacitly assume that dealing with elements in state space is 
0(1), which is strictly speaking only guaranteed to be true in a finite dimensional 
setting. However dealing, e.g., with curves on the real line numerically can still 
be of 0(1) if only the relevant parts of the curve are actually calculated. Hence 
the total complexity C of the method is d^nK, where d^ is a constant. Given an 
accuracy e the following inequality has to hold true additionally, 

di , d 2 . 

1 p= < e, 

n s y/K 

whence we end up with the simple constraint minimization problem to minimize 
the complexity 

C = d^nK — > min 

given the previous inequality on accuracy. Its asymptotic solution is given by 
C = Oie- 2 - 1 / 3 ) with K = 0(e" 2 ) and n = 0(e~ 1 / s ). This can be improved by 
multi-level methods [27J EH H] to a complexity estimate of order almost 0(e~ 2 ), 
which is in turn the complexity of one dimensional MC integration. In other words, 
the complexity is equal to the integration of a functional where evaluating at a 
single point is of order 0(1). Multi-level methods improve by telescoping errors on 
different levels of discretization n. However, in this case the asymptotic complexity 
is not improved by higher-order methods anymore, since it depends on weak and 
strong convergence orders so that one is restricted to low order Euler-like methods. 

If we perform the same complexity analysis in case of higher order discretization 
schemes with a QMC algorithm instead of an MC algorithm we obtain an asymp- 
totic complexity C = O^- 1 - 1 / 3 ) with K = 0(e _1 ) and n = 0{e^/ s ), which is 
indeed considerably better than multi-level MC in case of higher order methods 
(s > 1). On the other hand it is not better than multi- level QMC [TH] theoretically 
could be. We emphasize that a multi-level QMC is theoretically far from being 
understood, and additionally we would need strong order 1 methods which are not 
always at hand. An additional problematic aspect is the need of high-dimensional 
integration spaces, where QMC is not known to perform well anymore. 

We claim that standard QMC with high-order weak approximation schemes is 
superior to multilevel MC due to the low dimensionality of the integration space 
R d ™ as long as the QMC order of convergence is understood to hold true. For accu- 
racy e the dimension of integration space is of order 0(e _1 / s ), which in real world 
implementations is often sufficiently small such that the QMC order of convergence 
is ensured. 

The goal of this work is therefore twofold: first we want to show how actually 
the theory of weighted spaces applies to the HJM equation. We even show that we 
have a simple weak approximation method of order 2 within this setting. Second, 
we claim that a QMC algorithm integrating the resulting functional is numerically 
efficient. We underline this statement by a calibration of a time-homogeneous, 
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non-linear, diffusive HJM-equation to caplet data, i.e., we calibrate this equation 
to ten volatility smiles simultaneously. Our method is not only fast enough for the 
calibration of the model, but also the computer programming itself is almost as 
easy as a standard Euler-Maruyama scheme due to the use of a splitting approach. 

Let us compare our results to well-known and recent results on splitting schemes 
and weak approximation methods for SPDEs. In contrast to classical results on 
the Lie-Trotter splitting such as [H El HZ1 E IH EH 133 H21 I2H1 [201 HI], we 
focus on a higher order method for nonlinear problems in the spirit of [38] . hence 
allowing us to conclude the practical efficiency of the method as explained above. 
The topic of weak approximation for SPDEs was recently analysed in [9] . While 
there, the focus was on space-time white noise driving the system, we consider 
only finite-dimensional noise, and can obtain the same rate of convergence as in 
the finite-dimensional setting with bounded and smooth vector fields. Contrary to 
[36] , our model is inherently infinite-dimensional and does not allow a reduction to 
a low-dimensional stochastic differential equation. 

2. Weighted spaces and analysis of stochastic partial differential 

equations 

We provide an overview of the theory of weighted spaces that is at the core of 
the presented numerical method. For more details, see also [40l [12l [TT1 ITO] . 

2.1. The generalised Feller condition. Given a fixed t > 1, we consider the 
following setup. 

(1) For i = 0, • • • ,£, (Hi, (•, -)Hi) is a separable Hilbert space, and its norm is 
denoted by 1 1 x 1 1 j/^ := (x, x} 1 ^ 2 . 

(2) Hi + i is compactly and densely embedded into Hi for i = 0, ...,£ — 1. 

(3) A: dom ^4 C Hq —> Hq is the generator of a strongly continuous semigroup 
of contractions (S t )t>o on Hq. 

(4) For i = 0, ...,£- 1, A: H l+1 -> H t is bounded. 

(5) For i = 1, . . . ,£, (S t )t>o is a strongly continuous semigroup of contractions 
on Hi. 

In many cases, it will be adequate to choose Hi :— dom A 1 , e.g., if A is a differential 
operator on a bounded domain. If, however, A is a differential operator on an 
unbounded domain, dom A will usually not be compactly embedded into Hq. As 
we are interested in the HJM equation, where the underlying space variable varies 
in [0, oo), we consider the above, more general setup. 

Definition 1. Let i € {0, ...,£}. Given a left-continuous, increasing function 
p: [0, oo) — > (0,oo) with lim tl _ i . 00 p(u) = +oo, set ipi( x ) '■— p(\\ x \\Hi), we define the 
enveloping space 13^(7^) := {/ e C k (Hi): ||/||^ 4 ,fc < oo}, where C k (X) denotes 
the space of k times continuously Frechct diffcrcntiablc functions and 

k 

(1) ||/IU,fc := H/kj with l/l^- := sup M^W&fWhjm- 

Here, Lj(Hi) is the linear space of bounded multilinear forms a: — > R en- 
dowed with the norm 

(2) HUjCffi) := sup \a(xi,...,Xj)\, 

x-t,. --,Xj &Hi 
\\xi\\ Hi ,i=l,...,j 

which makes (Lj(Hi), \\-\\L-(Hi)) a Banach space. 
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Given an orthonormal basis (ej)j^, define the space A(H{) of bounded smooth 
cylindrical functions Hi — > M by 
(3) 

A(Hi) := {/: / = fl((-,ei) Hi ,.-- ,(-,e w )^) for some JV G N and g G C£° (K^) } . 

The closure of A(H 2 ) in Bf(Hi) is denoted by Bp (Hi), k > 0. 

Remark 2. The above assumptions on the weight function ^ are very restrictive. 
A weaker assumption on the weight function on which our analysis can be per- 
formed would be that the sets {x G Hi: ipi(x) < R} are weakly compact, and hence 
bounded, in Hi, and that ipi is bounded on bounded sets. This is applied in Sec- 
tion o 

Applying [?U1 Corollary 5.3, Remark 5.4], we see that our space Bp (Hi) coincides 
with the space WC^, i defined by M. Rockner and Z. Sobol. Hence, the following 
result is proved in 001 Theorem 5.1]. 

Proposition 3. There exists an isometric isomorphy from Bq (Hi)*, the dual space 
to Bq (Hi), to the space 

.~~ 

where the latter space is endowed with the norm H/iH^,* := J H . ipi(x)\fj,(dx)\, \fi\ 
denoting the total variation measure to \i. The inverse of this isometry is given by 
W) ■= I Si f(x)lM(dx) for all f G Bp(H t ) and p G M^(Hi). 

This result allows us to obtain a generalisation of the well-known Feller condition 
for the strong continuity of operator semigroups on Cq(D), D a locally compact 
topological space, to the infinite-dimensional setting. 

Corollary 4 (generalised Feller condition). Fix i G 0, . . . , I. Let (Pt)t>o be a family 
of continuous operators on Bp (Hi) satisfying the generalised Feller condition, i.e., 

(1) Po = I, the identity on Bp (Hi), 

(2) P t+S = P t P s for s,t>0, 

(3) \\Pt\\ Tl p,*i, n < C for all t G [0, e) with some C > and e > 0, where 

L(Bp (Hi)) is the space of bounded and linear operators on Bp (Hi) and 
endowed with the operator norm, and finally 

(4) lim t ^ + Ptf(x) = f(x) for all x G X and f G Bp(H t ). 

Then, (Pt)t>a is a strongly continuous semigroup on Bp(Hi), i.e., for every f G 
Bp(Hi), lim^ +||P t /-/|k=0. 

Proof. This is an easy consequence of Proposition [3] and [THl Theorem 1.5.8]: We 
only need to prove that lim t ^ + t(Ptf) = ?(f) for all / G Bp (Hi) and I G Bp (Hi)*. 
But by Proposition [21 1(f) = j x f(x)fi(dx) for all / G Bp (Hi) with some fj, G 
Ad^^Hi). As lirnt_>.o+ Ptf(x) = f(x) for all x G Hi, an application of Lebesgue's 
dominated convergence theorem yields the claim. □ 

Hence, in contrast to the weak continuity of Markov semigroups for infinite di- 
mensional stochastic equations [B] , the above result allows us to work with standard 
strongly continuous semigroups. 

Usually, the difficult part in verifying the generalised Feller condition for a given 
Markov semigroup (Pt)t>o is proving that P t (Bp(Hi)) C Bp (Hi). The following 
result can often be applied to this problem. 

Theorem 5. For k > and i = 1, . . . ,£, C^(H^i) C Bf(H t ) is dense. 
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Proof. Apply Proposition 1161 to obtain an orthonormal basis (e„) ne N of Hi-\ that 
is simultaneously orthogonal in Hi. Defining A(Hi) using (e n /||e n ||jj i )neN! we see 
that every / = g((-, ei/\\ei\\ Hz ) Hz , . . . , (■, e N /\\e N \\ Hi ) Hi ) £ A{Hi) can be extended 
to a smooth cylindrical function on as 

(5) 5((-,ei/||ei||fl- 4 )j? 4) . ■ • , (■, e^/||ejv||fl-<)fl-<) 

= fl'(||ei||ir 4 ( , > e i)^i-i)-- • > ll e w||H s (-,eAr) i f ! _ 1 ). 

Whence A(Hi) C C£ (#<_].). 

Next, we show C£(i^_i) C Bf(Hi). Given / £ C£(Fi_i) and e > 0, we 
shall construct f e £ A(Hi) such that ||/ — /ell^fc < e. Let ttn denote the ifj-i- 
orthogonal projection onto spanje^ : j = 1, . . . , N}. For i? > arbitrary, we esti- 
mate 

k 

||/-/°7rjv||^<Xl SU P ^(*)~ 1 H^/(*)--D , 7(^)ll^(H0 
,-=i „ x „ eH - 
\W\\ Hi <R 

(6) +^)- 1 ||/-/°7rJv||cJ(i ?i) 

As 7Ttv : flj — > -Hi— l is of operator norm one for all N £ N, it is easy to see by the 
properties of p that the final term goes to zero as R goes to infinity, and hence can 
be made smaller than e/3 by choosing i? e , depending on / but not on N, large 
enough. For the first term, note that B :— {x £ H^. \\x\\ff t < -Re} is precompact in 
Hi-i. Hence, there exists 8 > such that \\D^ f{x) — f(y)\\ Lj ^ H .- j < e/3 for x, 
y £ B with \\x — yWn,-! < 8. Choose Ng according to Corollary IT71 to obtain that 
\\x — ttnxWhi-! < 8 whenever N > Ng and x £ B. 

Finally, choose f e : H l , Ns -> R in such a way that X^=o su PxeH, Ns \\D j f s (x) - 

M\ Hl 

D 3 (foir Ne )(x)\\ Lj(Hi) < e/3 and ||/ E || c *(j? iiJV5 ) < ll/ ^llc»(J?«,^)- Here > H i,N, ■= 
span {ej : j = 1, . . . , Ng}. Such a choice is always possible, as Hijy s is finite dimen- 
sional and we can thus apply a standard mollifying argument. It follows similarly 
as above that ||/ e — fo ttn s \\^ it k < and plugging the results together, we obtain 

(7) ||/-MU,*<e. 

Thus, A(Hi) C CliHi-x) C Bf ' {Hi), and the claim follows. □ 

Theorem 6. Fix i £ {1, ...,£}. Let (x(t, xo))t>o be a time homogeneous Markov 
property on the stochastic basis (fi, J 7 , (_? r t )t>o, P) with values in Assume that 

(1) the mapping — > Xo i-> a;(t,Xo) is almost surely continuous with 
respect to the norm topology on Hi for every t > 0, 

(2) i/iEo € flj and £ > 0, then x(t,xo) £ Hi almost surely, 

(3) for some e > and C > 0, E[^i(a;(t, Xo))] < Cip(xo) for all xq £ Hi and 
t £ [0, e], and 

(4) (x(t, xo))t>o has almost surely cadldg paths in the weak topology of Hi. 

Then, P t £ L(Bt{X)) for all t > 0, where P t f(x ) := E[/(aj(t,a:o))], (P*)t>o 
satisfies the generalised Feller condition, and hence, (Pt)t>o is a strongly continuous 
semigroup on Bq(X). 
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Proof. First, we prove that limt_>.o+ Ptf(xo) — f{xo) for fixed / e £>g l (i^) and 
x e Hi. Let R > \\x \\ Hi . Set B R := {x e H, : \\x\\ H% < R}, then 

\P t f(x ) - f(x )\ < E[\f(x(t, x )) - f(x )\] 

<E[\f(x(t,x Q )) ~ f(x )\xB R (x(t,x ))} 
+ E[\f(x(t,x ))\x Hi \B R {x(t,x ))} 

(8) +\f(x )\¥[\\x(t,x )\\ Hi >R] 

with xa{x) := 1, x € A, otherwise the indicator function of the set A. The 
Markov inequality yields 

(9) P[\\x{t,x )\\ Hi >R]< plR^Rtyitxfrxo))], 
and this term goes to zero as R goes to infinity. Furthermore, 

(10) E[\f(x(t,x ))\ XHi \B R (x(t,xo))} < \\fU ifi nMx{t,xo))x Hi \B R (x{t,x ))}, 

and dominated convergence proves that this also goes to zero as R goes to infinity. 
Finally, note that is weakly continuous, as / 6 £>Q ! (iii). Weak compactness 

of Br yields that \f{x) — f{xo)\ < 2 sup^g^^ < oo for x <G Br, and monotone 

convergence proves lim t ^o+ xq)) — f{xo)\xB R (x{t, xq))] = 0. Hence, we 

have shown lim t ^ + Ptf(%o) — f(xo)- 

Next, note that Pt(Cb(Hi^i)) C Cb(Hi-i), which is a consequence of the as- 
sumption of almost sure continuity of the mapping xq h > x(t, xq). As 

||ft/lk,o= sup i)i{x )-^\J{x{t,x a ))\\ 

(11) < ||/IU,o sup ^{xor^iixit, x ))} < CH/11^,0 

for t G [0,e], this proves by Theorem [3] that P t e L(Bp(Hij) for < G [0,e]. As 
the semigroup property is satisfied due to the Markov property, an induction shows 
P t € L(Bp(Hi)) for all t > 0. Thus, Corollary H proves the claim. □ 

2.2. Application to stochastic partial differential equations. Let (x(t, £o))t>o 
be the solution of the stochastic partial differential equation 

d 

(12a) dx(t,x ) = (Ax(t,x ) + V {x(t,x )))dt + ^rVj{x(t,x )) o dWl , 

3 = 1 

(12b) x{0,x o ) = x o . 

Here, (Wt)t>o is a d-dimensional standard Brownian motion. The vector fields Vj 
are assumed to be of the form Vj(x) — gj(Lx), where gj £ C^°(R Ar ; Hi) is a smooth 
function on with values in Hi, and L: Hq — > M. N is a bounded linear mapping. 
These are typical assumptions for HJM models to be applied in practice, see [15] . 
Then, it follows that (IT21 admits unique solutions in every space Hi, i = 0, . . . ,£, 
given that the initial value xo is smooth enough. 

Lemma 7. Fix (5 > and i € {0, ...,£}. For some e > 0, there exists C > such 
that 

(13) E[cosh(/3||a;(i,a;o)||ff 4 )] < C cos\i(P\\x \\h,) for xq e Hi andte [0,e]. 
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Proof. We apply Ito's formula. For m > 2, 

d\\x(t, a; )|||^ = m\\x(t,x )\\%™~ 1 '{x(t,x ),<ix(t,Xo))iii 

+ ^m((m- l)\\x(t 7 x )\\h™~ 2) (x(t, x ), dx(t, x )) 2 H . 

+ \\x(t, xo)\\ 2 ^ 7 j n ^ 1 ' > (dx(t, x ),dx(t,x )) Hi ^ 
= m\\x(t, xo)||# 8 TO_1) ({x(t, x ),Ax(t, x ))ff s dt 

+ (x(t,x ),Vo(x(t,xo))} Hi dt 

d 

+ ^2(x(t, xo),Vj{x{t, x ))) Hi dWi 



1 d 

+ -m||.T(t,xo)|!^ m - 2) ^((m-l)(x(^.T ),^(x(t^o)))^ 



(14) 



\ x {% x o)\\h z ( v j( x {t, x o)), Vj ( x (t, x ))) Ht j dt 



Taking expectations, the boundedness of the Vj and the dissipativity of A yield, 
as all moments are uniformly bounded by Theorem 7.3.5], a constant C > 
independent of m > 2 such that 
(15) 

E[||x(t,x )||l£ l ] < \\x \\% + Cm [ E\\\x( S ,x )\\ 2 ^- 1 +m\\x(s,x )\\^ m - 1) 

Jo 1 

For m = 1, we similarly obtain 

(16) E[||z(t,z )|iy < \\x \\%. +C [ E[||x(s,x )|k + l]ds, 



ds. 



and trivially, E[||x(i, xo)||°] = 1. Note that cosh(w) = Ylm=o (2m)! • Summing up, 
the monotone convergence theorem proves 

E[cosh(/3||x(i,xo)lk)] < cosh^Haolk) 



+ Cf3 / E 



2m-l 



m=l V 



+/ 3 E7Sv/ 32(m " 1) n-(^o)ii^r 1) 

m=l 1 j - 



<cosh(/3||x || ff J 



/3 



(17) 



E[sinh(/?|jx(s,xo)|| ff! ) 
/3cosh(/3||x(s,x )||j? i )]ds. 



< 1 for m > 1. 



Here, we have used that sinh(u) = X)m=i (2m-i)l ' an< ^ that j m _ ] 
As sinh(u) < cosh(u), we obtain that with a constant C > depending on /3, 



E[cosh(/3||x(f, x )\\ Hi )] < cosh^llxollnj 



(18) 



C 



1 [cosh(/3||x(s,x )||H i )]ds. 



The method of the moving frame (see [13]) allows us to conclude that 
E[cosh(/3||x(i, xo)||/fJ] < oo for t > 0. Hence, an application of Gronwall's inequal- 
ity proves the claim. □ 
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Hence, the choice of weight function ipi^(x) :— cosh(/3||a:||// i ), j3 > 0, is appro- 
priate. This is particularly important in the application of our results to the HJM 
equation, see Section [5] 

Corollary 8. Given i g {1, ...,£} and f3 > 0, the Markov semigroup (Pt)t>o of 
(x(t, Xo))t>o is strongly continuous on B '' fi (Hi). 

Proof. Under the given assumptions, we can prove the conditions of Theorem [6] 
using [I Theorem 7.3.5]. □ 

Choose some £o <E {1, . . . , £} and (3q > 0. We perform an analysis of the infinites- 
imal generator Q with domain domG of (Pt)t>o, considered as strongly continuous 
semigroup on B^ e o.f>o (Ht ). In the following, Vf(x) := Df(x)(V(x)) denotes the 
directional derivative for sufficiently smooth functions / : H, t — > R and vector fields 
V-.Hi^Hi. 

Lemma 9. Fix i 6 {0, ...,£}. For j — 0, . . . , <i and f G A(H t ), Vjf € A{Hi). 
Furthermore, the directional derivative f i-> Vjf defines a bounded linear operator 

from tit'* (Hi) to B^f (fli), k > 1. 

Proof. The special form of Vj proves Vjf E A(Hi) for / € A(Hi). The estimate 
||V^/||^ i(3 < C||/||^, i/3 can be shown by a direct calculation using the boundedness of 

Vj and its derivatives, and the result follows from the density of A(Hi) in B^' 13 (Hi). 

□ 

Lemma 10. Fix i <E {1, ■■■,£} and fi\ < /3 2 - The operator f i-> Df(-)(A-) maps 
A(Hi-i) to A(Hi), and defines a bounded linear operator from B^~ 1,fl1 (Hi-i) to 
Bi^(H % ),k>l. 

Proof. Given / 6 A(Hi-\), there exists a iZi_i-orthogonal projection 7r with finite- 
dimensional range such that f on = f. Hence, Df(x)(Ax) = Df(x)(nAx), and it 
is easy to see that this function is in A(Hi). The boundedness is again shown by a 
direct calculation, where we apply that u cosh(/3iw) < C cosh(/3 2 w) for all u € [0, oo) 
with some constant C > 0. □ 

An application of Ito's formula, see [HI Theorem 7.2.1], yields that for i € 
{!,...,£} and /e^(ffi-i), 

d 

(19) 0/(a!) = fl/(i)(^) + (7 /)(i) + -X;^/)(i) forxGHi. 

i=i 

Theorem 11. Fiic i 6 {1,...,^}. For j > and < /?i < /? 2 , fie operator 
Q: B^^Hi^i) — > B^ P2 (Hi), given by the right hand side of (TT9"]) . is well-defined as 

a bounded linear operator. Furthermore, for f3 £ (0,/3q), B 2 '°~ 1 '' 3 (H(, -x) C domC?, 
and on i/iis space, Q = Q. 

Proof. The boundedness of follows from Lemmas [5] and [TO] For the second 
property, note that Q maps B 2 e °~ 1,fl (He. -\) into B^"' 13 " (Hi ) as a bounded linear 
operator, £/ = <5 on ^l(i?£ _i), and that Q is a closed operator. Hence, Q = Q 
follows from a density argument. □ 

Corollary 12. Fix (3 € (0,/3o). Given k e {0, . . . — 1}> we ^ aue Taylor 
expansion 

(20) p t / = jflf + R ^f f° r f e < fe °;i) +1> '" (^o-ck+i)). 
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where the operator Rt.k'- ^2(fc+i) +1> C^?a-(fc+l)) ~^ B^ e °> p ° (H(. ) is bounded uni- 
formly in t 6 [0, e] /or given e > 0. 

Proo/. Theorem E] proves that Bj^g" 1 " '" (i^ _ (fc+1) ) -> (H lo ) is a 

bounded linear operator for J = 0, . . . , fc + 1. Hence, a standard Taylor expansion 
argument can be applied to prove the stated theorem. □ 

Lemma 13. For j3 > 0, k > and i e {1, . . .,£}, P t : Bf* (Hi) -> B^^ (Hi) is 
a bounded linear operator. Its operator norm is bounded uniformly for t £ [0, T], 
where T > can be chosen arbitrarily. 

Proof. This is consequence of smooth dependence on the initial value in 

By considering the sensitivity equations, see [SJ Theorem 7.3.6], all derivatives 

D 3 XQ x(t, xo)(hi, . . . , hj) are shown to satisfy bounds of the type 

(21) E[\\Di x(t, xo)(hi, . . . , hj)\\ p Ht ] < C p (||/n|| ff4 • • • WhjW^f for p > 2, 

where C' p is independent of xq. The boundedness of Pt in the norms given above 
then follows from the Cauchy-Schwarz inequality together with the property 
ccosh(2u) < cosh(u) 2 < Ccosh(2u) for some constants c, C > 0. Due to Theo- 
rem HH it follows that P t (A(Hi-t)) C B^!' ,l3 (Hi). A density argument proves the 
claim. □ 

Corollary 14. For k > 2, i e {0, . . . Jo - 1} and f3 e (0, f3 ), B^{Hi) is a core 

for g. 

Proof. Applying (T3l Proposition II. 1.7] , this is clear from Lemma 1131 as B^ 1,13 (Hi) C 
domQ is invariant with respect to the semigroup and dense in Bq'^ (Hi). □ 

3. The rate of convergence of splitting schemes for stochastic 
partial differential equations 



As numerical discretisation scheme, we suggest the use of a splitting scheme. 
Decomposing the drift coefficient further, Vq — 2 TO= iVo,mj we define the split 
problems 

(22a) ^-x ,o(t, x ) = Ax ,o(t, x ), 

at 

(22b) -r:Xo,m(t,xo) = Vo, m (xo, m (t,x )), m=l,...,M, 

at 

(22c) dx j (t,x ) = V j (x j (t,x ))odWi, j = l,...,d. 

We stress that all of these problems can be solved by finding the corresponding 
deterministic flows; in the case of j = 1, . . . , d; we need to evaluate the flow induced 
by the vector field Vj at the stochastic time Wl . In particular, the processes 
(xo,m(t, Xo))t>o, m = 0, . . . , M, are deterministic. The split semigroups are defined 
by P?' m f(x ) := f(x , m {t,x )) and P] f(x ) := E[f( Xj (t,x ))], j = l,...,d. We 
consider the following splitting schemes. 

Lie- Trotter splitting, forward ordering: The Lie- Trotter splitting with 
forward ordering is of first order and reads 

(23) Q\l^ d f := P°<° P°J . . . P°f Pl t . . . Pi t f for / G B^o { H £o ). 

Lie- Trotter splitting, backward ordering: The Lie- Trotter splitting with 
backward ordering is obtained by reversing the order of the operators in 
the Lie- Trotter splitting with forward ordering, 

(24) Q(At) Wd / := Pit- PltP°At M • • • P°JP Atf for / e B^o (fl^), 
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and is also of first order. 
Ninomiya-Victoir splitting: The Ninomiya-Victoir splitting is a generali- 
sation of the well-known Strang splitting to more than two generators and 
reads 

n NV f ._ 1 p 0,0 / p 0,l p 0,M p l prf 
^(At)J ■— 2 At/2 I •'At ■ • • r A( 'At ■ • • 'At 

(25) + Pl, . . . PiA M . . . P°j)P° A ? /2 f for / £ (H ia ). 
It is of second order. 

The theory of Section [2] now applies not only to the continuous semigroup (-Pt)t>o> 
but also to every split semigroup (P t ' m )t>o and (P/)t>o, yielding spaces invariant 
to the dynamics of (Pt)t>o on which we can apply the generators £/, £?o,m an( i Sj, 
m = 0, . . . , M and j — 1, . . . , d, and 

M d 

(26) = X^°> m + X^- 

m— i — 1 

Hence, we obtain the following result. 

Theorem 15. Let (3 £ (0, /?o); and assume that (Q(At))At>o ts an V splitting 
approximation of (Pt)*>o ^aserf on i/ie spZii semigroups (P* ' m )t>o anc ^ (P/)t>o> 
m = 0,...,M, j = 1, which is of formal order s G {1, . . . ,£o — 1}- P° r 

/e^)(ffo), 

(27) ||P t / - Q£ /n) /|k +1 ,, < C T n- 5 ||/||^ „ 2(s+1) . 

Proof. The theory of [35] yields this result in the following manner. Clearly, all split 
semigroups are stable on the space Bq B+1,f3 ° (iJ s+1 ) in the sense that the operator 
norms of the operators are bounded by exp(Ci) with some constant C > indepen- 
dent of t and of the semigroup. Furthermore, for every t > 0, Pt is a bounded linear 
operator on fifu+i)^ ) ^ v Lemma [TBI and on this space, we have that all gener- 
ators of the split semigroups and the original semigroup are well-defined together 
with their products, and satisfy 

(M d \ a 

m=0 3=1 ) 

Hence, we obtain the claimed result from [23} Theorem 2.3, Sections 4.1, 4.4]. □ 

4. Symmetrically weighted sequential splitting 

Applying the theory of [23J [23] allows us to obtain asymptotic expansions for 
the forward and backward ordering of the Lie- Trotter splitting and the Ninomiya- 
Victoir splitting if the function / is sufficiently smooth. Using symmetry, we can 
even prove that the Ninomiya-Victoir splitting and the symmetrically weighted se- 
quential splitting, going back at least to |41[ equation (25)] and given by 

(n n \ nSWSS f 1 / / n LTfwd\n f , / n LTbwd\n A 

Wt,n J ■— 2 \W(t/n) ) J + W((/n) ) J ) i 

have asymptotic expansions not only in n _1 , but even n~ 2 . Hence, every extrapo- 
lation step would improve convergence by two orders. In particular, the symmet- 
rically weighted sequential splitting is of second order. Comparing the dimension 
of integration space of different second order schemes and in view of possible ex- 
trapolations we use SWSS in our numerical computations detailed below. Indeed 



EFFICIENT NUMERICS FOR HJM MODELS 



11 



dimension of integration space for the Ninomiya-Victoir scheme is n(d+ 1), whereas 
sequential splitting leads to dimension nd + 1 . 

5. Application: the Heath- Jarrow-Morton equation 

As application of our theoretical results, we provide a numerical method for the 
efficient simulation of the Heath- Jarrow-Morton equation of interest rate theory, 
ft is of the form specified in (|12p . where the infinitesimal generator is given by 



the differential operator xr. In order to include a stochastic volatility process, the 
Hilbert spaces Hi, i = 0, . . . ,£ are specified as follows. We set 

Hi := {/i € L 1 1 oc ((0, oo)) : h is i + 1 times weakly differentiable and 

(30) h',. . . M t+1) € l4((0,oo))} x R. 
Here, < ceo < ■ ■ ■ < a i , and 

(31) l£((0,oo)) := \ heLl c ((0,oo)): f r{x) 2 exp(ax)dx < oo\ . 

{ -Ao,oo) J 

It is easy to see that Hj+i C Hi for i e {0, ...,£— 1}, and that every function in 
is continuous and bounded (see also [II]). The scalar product on Hi reads 

((/ii,wi),(fc2,V2))fli := ^i(0)ft 2 (0) 

(32) +^ / /i^ m '(a;)/i2"^(a;) exp(aix)da; + -yiu 2 . 
With the induced norm, 

(33) A-.H^^H,, {h,v)^ (h',-av), 

where a > is a constant, becomes a bounded linear operator. It agrees with the 
generator of the shift semigroup on the first component of Hi, i £ {0, . . . ,£ — 1}. 

Consider the Heath- Jarrow-Morton equation with stochastic volatility in Ito 
form, 

dr(t,r Q ,v Q ) = (Ar(t,r ,v ) + auju(r(t,r ,Vo),v(t,v )))dt 

d 

(34a) + J2 °j (<t, r , «o), v(t, u ))dW/ , 

i=i 

(34b) d«(t,« ) = -av(t,v )dt + ^TjjdWi , 

3=1 

(34c) r(0, r ,w ) = r , 

(34d) »(0,«o)=«d. 

The stochastic volatility uq) is chosen as a mean-reverting Ornstein-Uhlenbeck 
process. The HJM drift satisfies the condition 

(35) QiHJM (h, v) (x) = y~] aj (h, v) (x) I CTj-(/i,w)(£)d£. 

We assume that aj are of the form required in Section l2~2l i.e., Oj (h, v) — Qj(Lh, v), 
where g € Cf(R N+1 ; H e ) and L: H R N is bounded linear. Rewriting the 
equation in Stratonovich form, we see that 



1 d 

(36) Vo(h,t>) = anjM(h,v) - ^^Daj(h, v){<Jj{h, v)), 



3=1 
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and it follows easily that Vo(h,v) — go(Lh,v) with some go £ C^°(R Ar+1 ; Hi). 
Hence, Theorem [15] applies to prove the optimal rate of convergence of s of a 
splitting scheme for sufficiently smooth functions / : Hq — > R, given that the initial 
value satisfies (ro,wo) £ H s+ \. 

5.1. The money market account. In order to calculate standard payoffs, we 
not only need the instantaneous forward curve, but also the money market account 
(B t )t> Q . It is given by B t = exp(z(t, 0)), where 

(37) dz(t,r Q ,v ,z ) = R t dt, z(0, r , v , z ) = z , 

and can therefore be easily included into our splitting scheme. Here, we denote by 
R t := r(t,ro,vo)(0) the short rate induced by our HJM model. 

To recover the optimal rate of convergence, we argue as follows. On the product 
space Hi := Hi x R, we consider the weight function i(>i t p(h, v, z) := ipi^(h, v) + z' 2 
(see Remark [2]). As proved before, 

(38) E[V>i,/9(r(t, r ,v ),v(t, v ))] < exp(Ct)ipi,p(r , v ). 
Furthermore, as Rt < ipi,p(T(t, ro, vq)), 

E[z(t,r ,v ,z ) 2 }<z 2 +f E[z(t,r ,v ,z ) 2 }ds+ [ E[R 2 t ]ds 



o 



(39) <z Q + E[ipi t p(r(t,ro,vo),v(t,vo),z(t,ro,vo,z ))]d 



s. 



o 



Altogether, an application of Gronwall's inequality proves 

(40) E[ipi^(r(t, r ,v ), v(t, v ), z(t, r Q , v , z ))] < exp(Ct)ijji } p(r , v , z ), 

and we can apply the above theorems to all functions contained in (Hi x R) 
by evident modifications of the above proofs. 

Now, the money market account is not included in the above setting. More 
precisely, B t = exp(#(t, ro, Vo, 0)), and this growth is larger that the quadratic 
growth admitted by ijjijj. We deal with this problem in the following way: Actually, 
z(t, ro, vq, 0) should be nonnegative from an economic point of view. Hence, we 
replace the money market account by B t := cxp($(z(£, ro, vq, 0))), where $ : R — >• M 
is C°° with bounded derivatives, satisfies &(z) = z for all z > —K, and is bounded 
from below by —2K with some K > 0. In our numerical experiments, performed 
using the model calibrated to the data from |30| . we never encountered paths with 
z(t)0) < 0. Furthermore, even if z(t, 0) becomes slightly negative on some paths, 
this is numerically innocent, as we can adjust K accordingly. We want to stress that 
our modification only acts on economically dubious paths where the money market 
account falls significantly in the long run, and neither limits temporary decrease, 
or any increase whatsoever. 

Clearly, B^ 1 < exp(2K). Hence, the modified payoff of a zero coupon bond with 
time to maturity S, 

(41) f(h,v,z) :=exp(-$(z))exp(- f h(s)ds), 

Jo 

is included in our setup, and lies in B^ ' 13 (Hq x R) for all k > if (3 > is chosen 
large enough: first, note that / depends on h only via the bounded linear functional 
£: H -» R, 1(h) := f* h(s)ds. It follows that 

(42) \f(h,v,z)\ < exp(2K)exp(C\\h\\no)- 
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Choosing (3 > C, the claim is proved, as it is clear that we can approximate / by 
functions of the form (h,z) i-> exp(— $>(z))(p(l(h)) with ip € Cg°(R) in the norm of 

sf ^(iJ xR). 

While standard payoffs, such as caplets and swaptions, do not satisfy the smooth- 
ness assumptions required in our results, we can at least prove that they are con- 
tained in a space on which convergence - albeit without rates - is ensured. A similar 
argument as for the bond price can be used to prove that the modified payoffs of 
caplets, 

(43) f(h, v, z) = exp(-$(«))(Li(ft) - K)+, 

where L$(h) := | (exp ( /i(r)dr) - 1) is the LIBOR rate, and payer swaptions, 



f(h,v,z) = 
= exp(— $(z)) 



Vexp(- / /i(r)dr) exp / 
,i=i Jo V \ J (i-K 



h{r)dr - (1 + SK) 



+ 



are contained in B^°- i3 (Hq x R). Here, however, taking the positive part makes 
these functions nonsmooth. As the space B^ffs+i) x ^) °^ mnc tions on which a 
rate of convergence is proved is dense in B^ - 13 (Hq xK), we still obtain convergence. 

6. Numerics for the Heath- Jarrow-Morton equation 

We present the results of numerical computations for a Heath- Jarrow-Morton 
model. We do neither claim that the chosen HJM model is particularly well suited 
nor that the chosen calibration strategy is the best. We only want to demonstrate 
that a non-linear infinite-dimensional HJM model with stochastic volatility can be 
efficiently calibrated to market data with a satisfactory result. 

First, a numerical calibration to caplet prices is performed, afterwards, a payer 
swaption is priced using the calibrated model. In our numerics, space discretisation 
is performed using piecewise affine and continuous functions, where the mesh is 
aligned with the time mesh. Hence, the partial differential equation 

d d 
(44) — r fi(t,r )[x) = — r 0fi (t, r )(x), r ,o(0, r )(x) = r Q (x), 

is solved exactly by shifting tq. 

6.1. Calibration. We demonstrate the efficiency of the presented method by per- 
forming the calibration of a parametrised, time-homogeneous Heath- Jarrow-Morton 
model to the caplet volatility surface provided in [30]. Note that the bond prices 
given there are automatically reproduced in our model by choosing them as the 
initial value. 

We set d = 3, and specify <jj(h,v) = gj(h,v)Xj. Here, \j is assumed to be of the 
exponential-polynomial type [2], \j(x) = J2T=o a j,i xl ex P( — P x )- It is easy to see 
that under such assumptions, the regularity required in [HJ Section 5.2] is satisfied. 
In our experiments, we choose io = 2. 

There are several economically sound possibilities for choosing gj . Guided by the 
Cox-Ingersoll-Ross model, one could choose gj(h, v) = y/\vh(tj)\ with some tj > 0, 
where the absolute values are necessary as we cannot guarantee positive interest 
rates by this approach. This ansatz, however, is not contained in our general setup, 
as gj is not a smooth function of h. 

Instead, we assume gj{h,v) — tanh(cj exp(?;) J* 3 h(s)ds). This ensures that the 
volatilities are bounded and vanish if the benchmark yields J* 3 h(s)ds driving the 
equation go to zero. We discretize the HJM-equation by the symmetrically weigthed 
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Figure 1. Calibration of the tanh-type volatilities 



sequential splitting scheme, as described in Section [4] The calibration is performed 
by combining a custom-written genetic algorithm, searching for global minima, 
with the Levenberg-Marquardt implementation from |34) to optimise locally. The 
model caplet values are calculated numerically, using 12 time steps per year and 
2048 quasi-Monte Carlo paths, based on the direction vectors for Sobol' sequences 
of Joe and Kuo [29] . 

All in all, 13 parameters are used to match 120 prices, and total calibration 
time is 14.5 minutes running on 16 cores of a Primergy RX200 S6 spotting 4 Intel 
Xeon CPU X5650 processor, each of which provides 6 cores. The calculation of 120 
option prices takes about .5 seconds and therefore merits to be called efficient. 

We are able to match the market volatilites taken from [30] very well using 
the tanh-type volatilies. Only the error in the earlier time slices is significant, 
see Figure [T] This is typical for models without jumps. These are well known to 
misprice options close to maturity. This behaviour can also be connected to the 
short end of interest rates depending more on announcements by central banks than 
random fluctuations. 

With respect to the martingale property of traded assets, numerical calculations 
show that bond prices and LIBOR rates satisfy the expected value property to a 
very high precision already using 2048 quasi-Monte Carlo paths. 

6.2. Pricing. As an application, we price an at the money payer swaption with 
a time to maturity of T = 5 years, where the underlying swap pays out quarter 
annually over three years, i.e., at the times Tj = T + iS for i = 1, . . . , 12 and 
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6 = .25. A reference computation with 16384 paths and 120 time steps per year 
yields the value 0.0281579. Using 2048 paths and 12 time steps per year, as in the 
calibration, we obtain 0.028074. The relative error is thus approximately .003. As 
the calculation of the coarser approximation takes .25 seconds, we have established 
the efficiency of the suggested method. 

7. Conclusions 

We introduce an analytic setup for the analysis of weak approximation methods 
for stochastic partial differential equations. The Heath- Jarrow-Morton equation of 
interest theory is shown to be included in the class where this approach is applicable. 
Moreover, the set of admissible test functions contains important payoffs such as 
caplets and swaptions. We argue that higher-order weak approximation schemes 
can be used together with QMC algorithms to obtain an efficient pricing method, 
which is even superior to multi-level MC. The efficiency of our numerical method 
is proved by the calibration of the model to given caplet data. 

Appendix A. Functional analytic results 

Proposition 16. Let (A, (•, ■) x), (Y, (■, -)y) be separable Hilbert spaces with norms 
\\-\\x and \\-\\y such that Y is compactly and densely embedded into X . Then, there 
exists an orthonormal basis (e n )„ e N d Y of X that is simultaneously orthogonal in 
Y . Furthermore, lim^^ooHenlly 1 = 0. 

Proof. By the Riesz representation theorem, there exists a bounded operator k : X — > 
Y such that 

(45) (kx, u)y — {x, y)x for all x G X and y e Y '. 

With l: Y — » X the compact embedding, we set K :— lk. K is clearly compact 
and also symmetric, as 

(46) (Kxi,x 2 )x = (kxi,kx 2 )y = (%i,Kx 2 )x- 

Thus, there exists an orthonormal basis {e n ) n eN of X and a sequence (A„)„ e N C K 
decreasing monotonically to zero such that Ke n — A„e„, and we see that (e„) n6 N C 
Y. We obtain 

(47) (e n ,e m )y = \~ x {Ke n , e m ) Y = X~ l {e n ,e m ) x = A~ x <!> nim for n, m 6 N, 

1/2 

whence (e n ) ne N is orthogonal in Y and || e„ || y — X n , and the claim is proved. □ 

Corollary 17. Under the assumptions of Proposition \TS[ let 7Tjv denote the X- 
orthogonal projection onto Xn '■= span{e„ : n € N}. Then, 

(48) lim sup \\y - 7r N y\\x = 0- 

\\y\W<i 

Proof. By Parseval's identity, 

oo oo 

\\v - knvWx = (2/' e «)x= (y,Ke n ) Y 

n=N+l n=N+l 
oo 

(49) < sup A„ > (y, \ J e n )y < Ajv+i||y||y> 

n > N n=N+l 
— 1/2 

where we apply that ||e n ||y = A„ and that (e n )„ 6 N is orthogonal in Y. As 
(A„) rl gN decreases to zero, the claim follows. □ 
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